/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

qdd: A Quantum Computer C++ Library constructed using BDD's

Copyright (C) 1999-2007 David Greve <TheBeaNerd@gmail.com>

This library is free software; you can redistribute it and/or modify
it under the terms of the GNU Library General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.

This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Library General Public License for more details.

You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
USA.

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

#include "qdd.h"

/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

 qdd Support Functions

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

int chose() {

  // The random choice function

  int val = 1 + (int) (100.0 * rand() / (RAND_MAX + 1.0) );

  if (val > 50) {
    return(1);
  }
  else {
    return(0);
  }

}

int BIT (const Integer & a, unsigned int i) {

  // Tests the i'th bit of an Integer

  return(testbit(a,i));
}

Integer rand(int size) {

  // Produces a random Integer with (up to) size bits

  Integer r(0);

  int i;
  for (i=0;i<size;i++) {
    if (chose()) {
      setbit(r,i);
    }
  }

  return(r);

}

Integer bitpown(const Integer & a, unsigned int i, const Integer & n) {

  // Computes a^(2^i) mod n

  Integer r = a;

  unsigned int x;
  for (x=0;x<i;x++) {
    mod(sqr(r),n,r);
  }

  return(r);

}

Integer powmod(const Integer & a, const Integer & x, const Integer & n) {

  // Computes a^x mod n

  Integer r = 1;

  unsigned int maxBit = ((unsigned int) lg(x)) + 1;

  Integer p;

  unsigned int i;
  for (i=0;i<maxBit;i++) {
    if (BIT(x,i)) {
      mul(r,bitpown(a,i,n),p);
      mod(p,n,r) ;
    }
  }

  return(r);

}

Integer gcdm2(const Integer & a, const Integer & b) {

  // Thanks to John Cowles for the gcdm2 and gcdm1 functions
  // gcd(x,y) = x * gcdm1(x,y) + y * gcdm2(x,y)

  Integer Ak = a;
  Integer Bk = b;

  Integer Qk(0);
  Integer Akp1(0),Bkp1(0);
  Integer Ykm1(0),Yk(0),Ykp1(0);

  Ykm1 =  0 ;
  Yk   =  1 ;

  while (sign(Bk) != 0) {

    divide(Ak,Bk,Qk,Bkp1);

    Ak   = Bk;
    Bk   = Bkp1;

    Ykp1 = Ykm1 - (Qk * Yk);

    Ykm1 = Yk;
    Yk   = Ykp1;

  }

  return(Ykm1);

}

Integer gcdm1(const Integer & a, const Integer & b) {

  // Thanks to John Cowles for the gcdm2 and gcdm1 functions
  // gcd(a,n) = a * gcdm1(a,n) + n * gcdm2(a,n)

  // gcdm1(a,n) is used to compute the multiplicitive inverse 
  // mod n of a, where gcd(a,n) is assumed to equal (0 or) 1.

  Integer Ak = a;
  Integer Bk = b;

  Integer Qk(0);
  Integer Akp1(0),Bkp1(0);
  Integer Xkm1(0),Xk(0),Xkp1(0);

  Xkm1 =  1 ;
  Xk   =  0 ;

  while (sign(Bk) != 0) {

    divide(Ak,Bk,Qk,Bkp1);

    Ak   = Bk;
    Bk   = Bkp1;

    Xkp1 = Xkm1 - (Qk * Xk);

    Xkm1 = Xk;
    Xk   = Xkp1;

  }

  return(Xkm1);

}

Integer pmod(const Integer & n, const Integer & d) {

  // Comutes n mod d where "mod" is the "true" mod, not the C hack 
  // mod, % :)

  Integer r = n;

  qdd_assert(sign(d) > 0 , "pmod() : non-positive denominator");

  if (sign(r) < 0) {
    r = r + (1 - r/d)*d ;
  }
  if (r >= d) {
    r = r - (r/d)*d ;
  }

  return(r);

}

Integer invmod(const Integer & a, const Integer & b) {

  // multiplicitive inverse of a mod b.

  Integer r(0);

  qdd_assert(sign(b) > 0 , "invmod() : non-positive denominator");

  if (sign(a) == 0) return 0;

  return(((sign(b) == 0) ? 1 : pmod(gcdm1(a,b),b))) ;

}

void qdd_error(const char *msg) {
 
  // error .. die.

  cerr << "** qdd_error : " << msg << "\n" ;
  exit(1);

}

void qdd_assert(int test, char *msg) {

  // Assertion tester

  if (! test) {
    throw QuantumException(msg);
  }
}

int ilog2(Integer & arg) {

  Integer x = arg;

  // Number of bits required to represent x.

  qdd_assert( sign(x) > 0, "ilog2() argument less than 1");

  int i = -1;
  int p = 0;
  int n = 0;

  while (sign(x) > 0) {
    n |= p;
    p = BIT(x,0);
    x = x >> 1;
    i++;
  }

  return(i + n);

}

int ilog2(int x) {

  // Number of bits required to represent x.

  qdd_assert( x > 0, "ilog2() argument less than 1");

  int i = -1;
  int p = 0;
  int n = 0;

  while (x > 0) {
    n |= p;
    p = x & 1;
    x = x >> 1;
    i++;
  }

  return(i + n);

}

bdd Int2IndexBDD(int int_index, int max_index_bit) {

  // Converts an integer into a BDD of equivalent value.

  bdd bdd_index = bddtrue;

  int i;
  for (i=0;i<max_index_bit;i++) {
    if (int_index & 1) {
      bdd_index &= bdd_ithvar(i) ;
    }
    else {
      bdd_index &= bdd_nithvar(i) ;
    }
    int_index = int_index >> 1;
  }

  return(bdd_index);

}

bdd dag_satisfy(const bdd &M, const bdd &Fn) {

  // Constrains Fn to Fn' so that any variable choice from the
  // resulting function Fn' will satisfy the relationship M.

  bdd Mx = M;
  bdd Fx = Fn;

  while (! (Mx == bddtrue)) {
    
    qdd_assert((! (Mx == bddfalse)), "M cannot be satisfies");

    int x = bdd_var(bdd_support(Mx));

    bdd M1 = bdd_restrict(Mx, bdd_ithvar(x));
    bdd M0 = bdd_restrict(Mx, bdd_nithvar(x));

    bdd MCx = (! M1 ^ M0);

    Fx = bdd_ite(MCx,Fx,bdd_compose(Fx, ((! MCx) & M1), x));

    Mx = (!MCx) | MCx & M1;

  }

  return(Fx);

}

